To appear in The Astrophysical Journal 



Formation and Collapse of Nonaxisymmetric Protostellar Cores in Planar 
Magnetic Interstellar Clouds: Formulation of the Problem and Linear Analysis 

Glenn E. Ciolek 1,2 and Shantanu Basu 3 

vo > 

Q ; ABSTRACT 

o . 
<N ■ 

We formulate the problem of the formation and collapse of nonaxisymmetric pro- 
. tostellar cores in weakly ionized, self-gravitating, magnetic molecular clouds. In our 

formulation, molecular clouds are approximated as isothermal, thin (but with finite 
thickness) sheets. We present the governing dynamical equations for the multifluid 
system of neutral gas and ions, including ambipolar diffusion, and also a self-consistent 
' treatment of thermal pressure, gravitational, and magnetic (pressure and tension) forces. 

The dimensionless free parameters characterizing model clouds are discussed. The re- 

■ sponse of cloud models to linear perturbations is also examined, with particular empha- 
sis on length and time scales for the growth of gravitational instability in magnetically 

q ! subcritical and supercritical clouds. We investigate their dependence on a cloud's ini- 

tial mass-to-magnetic-flux ratio fiQ (normalized to the critical value for collapse), the 
dimensionless initial neutral-ion collision time f^o, and also the relative external pres- 
sure exerted on a model cloud P e xt- Among our results, we find that nearly-critical 
model clouds have significantly larger characteristic instability lengthscales than do 
more distinctly sub- or supercritical models. Another result is that the effect of a 

■ greater external pressure is to reduce the critical lengthscale for instability. Numerical 
simulations showing the evolution of model clouds during the linear regime of evolution 
are also presented, and compared to the results of the dispersion analysis. They are 
found to be in agreement with the dispersion results, and confirm the dependence of 
the characteristic length and time scales on parameters such as (jlq and P ex t- 
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formation 
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1. Introduction 
1.1. Background 

In recent years there has been a steady and significant accumulation of observational data that 
support the idea that magnetic fields play a pivotal role in the evolution of interstellar molecular 
clouds and star formation. For instance, polarization maps of extended regions in clouds reveal 
large-scale magnetic fields that in several instances are aligned with the short axis of clouds (e.g. 
Pereyra & Magalhaes 2004), which would be consistent with support from magnetic forces perpen- 
dicular to field lines. Infrared, far-infrared, and submillimeter maps of more localized regions show 
ordered field configurations with curvature, including "hourglass" shapes, that suggest dynamic 
interaction of magnetic fields and molecular cloud gas (Schleuning 1998; Schleuning et al. 2000; 
Fujiyoshi et al. 2001). 

Interferometric polarization maps of molecular cloud cores by Lai et al. (2001, 2002) indicate 
similar ordered field structures. Furthermore, using a modification of the Chandrasekhar-Fermi 
(CF) method (Chandrasekhar & Fermi 1953), they estimated magnetic field strengths and ener- 
gies that were large enough to exceed observed gas turbulence in their survey targets. A modified 
CF method was applied to the sub-mm polarization measurements of several prestellar cores by 
Crutcher et al. (2004), and they found that the mass-to-magnetic flux ratios of their cores were 
clustered about (within a factor ~ 2 above and below) the critical value for gravitational collapse. 
This result is consistent with earlier OH Zeeman determinations of the mass-to-flux ratios of proto- 
stellar cores (Crutcher 2004). Likewise, Curran et al. (2004) obtained sub-mm polarimetric images 
of two high-mass protostellar cores; a CF analysis of their data yielded mass-to-flux ratios for their 
two objects that were also nearly equal to the critical value for collapse. Finally, comparing CF 
results with OH Zeeman measurements for the immediate environment of the L184 prestellar core, 
Crutcher (2004) concludes that it may be an example of a nearly or barely critical core contained 
within a magnetically subcritical (i.e., supported) envelope. The grand picture painted by the vari- 
ous observational data above is one in which magnetic fields to a large degree control the structure 
and dynamics of star-forming molecular clouds, and has many elements of the original, unified 
theory of magnetically regulated star formation put forth early on by Mouschovias (1976, 1977, 
1978, 1979; see, also, the review by Mouschovias & Ciolek 1999). 

Ambipolar diffusion — the drift of neutral matter with respect to plasma and magnetic fields 
- initiates protostellar core formation in magnetically supported clouds. That is, because of im- 
perfect collisional coupling between neutrals and charged particles (including dust grains; Ciolek & 
Mouschovias 1993, 1994), gravitationally-driven diffusion of matter in the inner flux tubes of clouds 
redistributes mass and magnetic flux that leads to the formation and subsequent dynamical collapse 
of supercritical cores or fragments within massive, subcritical envelopes. This process was exhib- 
ited in the two-dimensional axisymmetric magnetic cloud simulations of Fiedler h Mouschovias 
(1993), and the axisymmetric thin-disk model calculations by Ciolek & Mouschovias (1994, 1995) 
and Basu Sz Mouschovias (1994, 1995a,b). Uniformly, these studies demonstrated the formation of 
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supercritical protostellar cores embedded in subcritical envelopes. They were also able to follow 
the resultant collapse of these cores well into the stage of dynamical infall, ending at the time when 
the central density was enhanced by a factor ~ 10 6 . 

Ciolek & Basu (2000; hereafter CB00) applied one of their axisymmetric magnetic disk models 
of core formation by ambipolar diffusion to the L1544 starless core. This model provided a good 
theoretical fit to the extended, subsonic infall profile as well as the distribution of mass in that core 
found by the earlier, complementary studies of Tafalla et al. (1998) and Williams et al. (1999). 
The predicted magnetic field strength in the CB00 model was later found to be in agreement with 
OH Zeeman observations of L1544 (Crutcher & Troland 2000; CB00). An analysis of submillimeter 
and millimeter maps and theoretical models of the dust temperature distribution in LI 544 were 
found to be consistent with the emission distribution calculated from the CB00 model (Zucconi, 
Walmsley, & Galli 2001). Subsequent detailed, multispectral line maps used to further examine the 
kinematics within the L1544 core have also shown reasonable agreement between measured infall 
speeds and the CB00 predictions (Caselli et al. 2002). Finally, a multispectral survey by Crapsi 
et al. (2004) has found that some of the infall features of the prestellar core L1521F (which, like 
L1544, is in the Taurus complex) are comparable to those of the CB00 model. 

However, despite the apparent successes of the aforementioned ambipolar diffusion models in 
describing the earliest stages of protostellar collapse, they are inherently limited because of their 
underlying assumption of axisymmetry. Axisymmetry is clearly a highly restrictive and idealized 
assumption, and is not likely to naturally occur in molecular clouds. In fact, sub-mm maps of star- 
forming clouds frequently indicate distinctively inhomogeneous and irregular large-scale structure 
and multiple cores (e.g., Andre, Motte, & Belloche 2001; Motte, Andre, & Neri 1998). Investigations 
at these wavelengths reveal irregularities on the scale of individual cores as well (e.g., Bacmann et 
al. 2000). Moreover, the observational studies of the L1544 core cited above find that it is definitely 
not axisymmetric. This is not a surprising result, as statistical analyses of data sets and catalogs 
of cores (e.g., Jones, Basu, &: Dubinski 2001; Jones & Basu 2002; Kerton et al. 2003) reveal that 
their shapes can generally be best fit with a distribution of triaxial ellipsoids, which are typically 
more oblate than prolate. 

Indebetouw & Zweibel (2000) simulated the formation of supercritical cores by ambipolar 
diffusion in infinitesimally thin, magnetic layers. Their models were two-dimensional and nonax- 
isymmetric, and they focused primarily on clouds that were initially subcritical. They followed the 
evolution of cores to when the maximum column density was a factor < 5 above the initial mean 
column density, and the irregular cores that developed were reminiscent of the sub-mm observa- 
tions cited above. 1 They also showed that the cores in their calculations developed on a timescale 



1 Some of the results in the interiors of their cores (scaling of magnetic field strength with central density, subsonic 
infall speeds, etc.) were similar to those found in the earlier axisymmetric model simulations over the same range 
of central density. This is also true for the nonaxisymmetric critical cloud model of Basu & Ciolek (2004 — see the 
denser, inner core regions in their Fig. 1). This agreement is probably the reason why the axisymmetric CB00 model 
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corresponding to the maximum linear growth rate of gravitational instability, and had a size equal 
to the corresponding wavelength (Zweibel 1998; see, also, § 3 below). Basu & Ciolek (2004; here- 
after BC04) presented numerical models displaying the gravitational collapse of multiple dense, 
asymmetric protostellar cores in thin (yet, not infinitesimally so) planar magnetic clouds. One 
of their model clouds had an initial reference state with a mass-to-magnetic flux ratio that was 
exactly the critical value (fio = 1, see §§ 2.2 and 3.1 below). Another model of theirs was initially 
supercritical by a factor of 2 (i.e., hq = 2). They found that, in the initially critical model, transfer 
of mass by ambipolar diffusion resulted in the eventual formation of supercritical cores contained 
within subcritical envelopes, similar to that which occurred in their earlier axisymmetric models 
that started with initially subcritical conditions. By contrast, the initially supercritical model of 
BC04 had a much more rapid and dynamical infall, with significantly greater infall speeds (sonic 
and supersonic) extending over much larger regions. Based on the physically distinct and different 
predictions of these two models, they noted that molecular clouds that are supercritical by more 
than a factor ~ 2 would be incompatible with the generally observed subsonic infall motions. 

1.2. Outline 

In this paper, we present the formulation for modeling the formation and nonaxisymmetric 
collapse of protostellar cores in planar magnetic clouds. In § 2 we describe the fundamental assump- 
tions and derive the necessary system of governing equations for a model cloud. The equations are 
put in nondimensional form, and the resulting free parameters of a model and their physical mean- 
ing are described. Our numerical method of solving the governing equations is also discussed. To 
provide a basis for understanding the underlying physics of ambipolar diffusion and gravitational 
collapse in magnetic sheets, the stability of model clouds is examined in § 3 by linearizing and 
Fourier-analyzing the governing equations. The results of the stability analysis are also compared 
with full numerical simulations of models in the limit of small-amplitude perturbations. In section 
§ 4 we summarize our results, and discuss their relevance to forthcoming fully nonlinear studies of 
nonaxisymmetric core formation and collapse. 

2. Physical Formulation 

We model clouds as isothermal thin sheets with temperature T, embedded in a hot and tenuous 
external medium of constant pressure P ex t- This simplifying assumption, the thin sheet approxima- 
tion, is based on the numerical simulations of ambipolar-diffusion-initiated formation of cores by 
Fiedler & Mouschovias (1993), who found that two-dimensional, axially symmetric magnetically 
supported model molecular clouds rapidly flatten and establish force-balance along the direction 
of magnetic field lines, even while contraction and dynamical collapse took place in the direction 
perpendicular to the field. These results were soon afterward incorporated in the protostellar core 



is a reasonable physical model for the interior of the asymmetric L1544 prestellar core. 
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formation studies of Ciolek & Mouschovias (1993; hereafter CM93) and Basu & Mouschovias (1994; 
hereafter BM94), who modeled molecular clouds as thin axisymmetric disks threaded by a vertical 
magnetic field, with hydrostatic equilibrium maintained along field lines at all times. Here we again 
adopt the formulation of CM93 and BM94, but now forego the restriction of axisymmetry. 

The appealing feature of the thin-sheet approximation is that it can be used to model the 
physical processes necessary to core and star formation in a theoretically tractable and realistic 
way, including a self-consistent treatment of both magnetic pressure and tension forces. Another 
advantage of thin-sheet models is that because the physical equations have been simplified by 
integration in one dimension (here, in the direction of the vertical magnetic field), the problem has 
been substantially reduced from its full three-dimensional complexity. Because of this, sheet models 
can provide significantly higher spatial resolution, and require much less computational time and 
storage than fully three-dimensional cloud models. There is also some observational justification for 
the use of the thin sheet approximation. The often observed alignment, or at least close correlation, 
between the projected magnetic field direction from polarization measurements and the core minor 
axis is evidence for flattening along the field (Basu 2000). Furthermore, an analysis of core shapes 
by Jones et al. (2001) and Jones & Basu (2002) implies that they are preferentially flattened along 
one direction. 

Of course, sheet models may not account for all of the observed morphological features of 
molecular clouds and their envelopes. Especially for those with large internal velocity dispersion, 
which can provide substantial support against self-gravity, and therefore be extended along magnetic 
field lines. It turns out though, that even in clouds or complexes with such substantial velocity 
dispersion or turbulence, thin-sheet models may still be used, so long as the application is restricted 
to dense sub-regions of the cloud. This has been demonstrated by Kudoh & Basu (2003, 2006), 
who analyzed the nonlinear support of stratified molecular clouds due to an ensemble of driven 
hydromagnetic waves. They found that because of density stratification, the largest, supersonic, 
velocities are in the low density envelope of a cloud, while a dense region near the midplane has 
transonic or subsonic motions. Hence, in this situation, thin-sheet models may reasonably be 
applied to a high-density fragment or sub-cloud, which will most likely be the site of subsequent 
star formation. 

The most significant systematic source of possible inaccuracy of a thin-sheet model is that it 
typically overestimates the strength of the gravitational field in the equatorial plane (i.e., the plane 
of the disk) when compared to less condensed or concentrated systems. For instance, the critical 
wavelength for gravitational instability in very thin non-magnetic clouds is found to be half the 
value of the critical length found for an equilibrium layer with an exponential atmosphere. This can 
be seen by comparing our equation (40) derived in § 3.1, in the limit of zero external pressure, to 
equations (13-36) and (13-42) of Spitzer (1978). However, the equilibrium layer calculation assumes 
an infinite vertical extent, which reduces the system's total gravitational binding energy; this is 
therefore an upper limit to the possible overestimate of calculating the gravitational field by using 
a thin disk. A better estimate was made by Kiguchi et al. (1987, see their Fig. 11), who showed 
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that the magnitude of the planar gravitational field for an infinitesimally thin disk exceeds that 
of a flattened cloud of finite extent by at most only ~ 40%, with the largest discrepancy between 
the two occurring for a very small region around the cloud center. BM94 presented a method of 
correcting for finite-thickness effects in the magnetic thin-disk approximation, deriving a quadratic 
correction (in terms of the local disk thickness) to the gravitational field. For the density regimes 
we intend to study (< 10 7 cm -3 ), the BM94 correction had the effect of altering the evolution of 
a model's core by about 10% during the early stages of core evolution, to < 30% at much later 
stages. Finally, we note that because the planar magnetic field in the sheet is calculated the same 
way that the planar gravitational field is (see eqs. [15] - [16b] and [20] - [21b] below), there is also 
a comparable overestimate in the magnitude of the magnetic tension force. Since the magnetic 
tension force opposes the gravitational force, the overestimate in the value of the gravitational 
field due to the thin-sheet approximation is offset by a corresponding overestimate of the planar 
magnetic field, thereby reducing the net effect of overestimating both fields. 

Thin-sheet models have been used by many workers in the field of molecular clouds and pro- 
tostars. They were used early on by Narita, Hayashi, & Miyama (1984) to study star formation 
in axisymmetric, non-magnetic clouds. As mentioned above, the thin-sheet approximation was 
developed and applied to ambipolar diffusion and core formation in isothermal magnetic molecular 
clouds by CM93 (who included the effects of dust) and BM94 (who included rotation and mag- 
netic braking). It was also later adopted by Li & Shu (1997a,b), who studied the equilibria and 
self-similar gravitational collapse of clouds with frozen-in magnetic flux (i.e., no diffusion). The 
approach to self-similar collapse during the later stage of core collapse with ambipolar diffusion 
in thin-disk clouds was described by Basu (1997). Ciolek & Konigl (1998) incorporated the ef- 
fect of the formation of a central gravitating point mass (a protostar) in numerical simulations of 
ambipolar-diffusion-driven dynamical core collapse in axisymmetric thin-disk clouds; a self-similar 
solution to this same problem was also provided by Contopoulos, Ciolek, & Konigl (1998). Tassis 
& Mouschovias (2005a,b) also investigated this particular topic, further extending the analysis by 
including a detailed calculation of multi-fluid effects on the conductivity of the infalling gas during 
the later stages of supercritical core collapse and accretion. 

A nonaxisymmetric 'toy model' idealization of a magnetically critical, flux-frozen turbulent 
sheet-like molecular cloud was suggested by Allen & Shu (2000). As noted in §1.1, Indebetouw & 
Zweibel (2000) first presented nonaxisymmetric simulations of core formation by ambipolar diffusion 
in infinitesimally thin magnetic sheets. BC04 presented models of nonaxisymmetric, gravitationally 
collapsing cores in magnetically critical and supercritical finite-thickness sheet-like molecular clouds, 
including ambipolar diffusion and its consequent effect on the dynamical evolution of clouds and 
cores. Li &: Nakamura (2004) and Nakamura & Li (2005) used the thin-sheet approximation to 
examine the combined effect of turbulent initial conditions and ambipolar diffusion in forming 
supercritical cores. 
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2.1. Fundamental Equations 

We present the necessary system of equations to model core formation in weakly ionized, 
magnetic interstellar clouds. As stated above, our model clouds are taken to be thin, with local 
vertical half-thickness Z(x, y, t) in a Cartesian coordinate system. By a sheet being thin we mean 
that for any physical quantity f(x,y,z,t), the criterion f/\V p f\ S> Z is always satisfied, where 
V p = xd/dx + yd/dy is the planar gradient operator. The magnetic field threading a cloud is 
taken to have the form 

B(x,y,z,t) = B ZtCq (x,y,t)z for \z\ < Z(x, y, t), (la) 

B(x,y,z,t) = B z (x,y,z,t)z + B x (x,y,z,t)x + B y (x,y,z,t)y for \z\ > Z(x,y,t), (lb) 

where -B 2i eq is the magnetic field strength in the equatorial plane of the cloud. In the limit \z\ — > oo, 
B — > B TC {Z, where B re f is a constant, uniform reference magnetic field very far away from the sheet. 
From now on, all physical quantities are understood to be a function of time t. 

Note that the condition on the planar gradient of physical quantities within the thin-sheet 
approximation implies that Z is the lower limit to the scales that can be described by our model. 
(For the gravitationally unstable modes, this condition is always satisfied, as the critical lengths 
always exceed this value — see eq. [40].) 

The unit normal vector to the upper and lower surfaces of the sheet is 

= ±^m,dx )x + {d Z, d y)y] 

[i + {dz/dxf + (dz/d y ) 2 \ 1/2 

Use of the integral form of Gauss's law yields the continuity equation for the normal component of 
the magnetic field across the upper and lower surfaces of the sheet, 

dZ dZ 
B z (x,y,±Z) - B x (x,y,±Z)— - B y (x,y,±Z)— = B z>eq (x, y). (3) 

The system of multifluid equations necessary to model a weakly- ionized, magnetic, self-gravitating, 
isothermal molecular cloud is given by equations (9a)- (9m) of CM93. (In this study we ignore the 
dynamical effect of interstellar dust grains. Hence, those terms representing grain contributions 
appearing in the system of fluid equations in CM93 are neglected.) As was done in CM93 and 
BM94, we simplify these basic equations in the thin-sheet approximation by vertically integrating 
them from z = —Z(x,y) to z = +Z(x,y). Doing so for the equation of mass continuity yields 

dcr n d d 

-^ + ^K^) + ^K^) = o , (4) 

where a Q (x,y) = p n (x,y)dz is the mass column density, and v n ^ x and v n>y are respectively the 
x- and y-components of the neutral velocity. In deriving equation (4) we have used the chain rule 
to obtain the velocity of the surface of the disk, 

dZ dZ dZ 

-dt =V ^ + V ^ ■ (5) 
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We have also employed a "one-zone" approximation, where we assume z-independence of physical 
quantities such as p n , and the planar velocity components v QtX , and v Qty . This approximation is 
also used for the planar components of the gravitational acceleration g x and g y that appear in the 
equations below. 

To derive the x- and y-components of the force equation (per unit area) for the neutrals, we 
use equation (2) and the total (thermal plus Maxwell) stress tensor 
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where 1 is the identity tensor, C = (k^T/rUn) 1 / 2 is the isothermal speed of sound; &b is the 
Boltzmann constant, and m n is the mean mass of a neutral particle (= 2.33 a.m.u. for an H2 gas 
with a 10% He abundance by number). Using equation (3), the symmetry conditions on B Zfiq and 
B z and the antisymmetry conditions on B x and B y about the equatorial plane, along with the 
divergence theorem, the integrated force equations are 
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(9) 



(-Fext + § GfJn) 

is the local effective sound speed. G is the gravitational constant. The second term on the right 
hand sides of equations (8a) and (8b) represent small modifications to the magnetic force due to 
planar gradients of the half-thickness Z. 

The equations for v hX and Vi tV , the x- and y-components of the ion velocity, are similarly 
obtained by vertical integration of the force equation for the ions. They are, respectively, 



1 T ni p 

Vn,x ~l "ma,g,x > 

+ — F 

O n 



(10a) 
(10b) 
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where the neutral-ion collision (momentum-exchange) time 

/ mi + m„ \ i 

r ni = 1.4 ^ . r • (11) 

V mi J ni(<7iy)iH 2 

The quantity m\ is the ion mass, which we take to be 25 a.m.u., the mass of the typical atomic 
(Na + , Mg + ) and molecular (HCO + ) ion species in clouds; (ow)\y{ 2 is the neutral-ion collision rate, 
and is equal to 1.69 x 1CP 9 cm 3 s _1 for H2 — HCO + collisions (McDaniel & Mason 1973). The factor 
of 1.4 in equation (11) accounts for the fact that the inertia of helium is neglected in calculating the 
slowing-down time of the neutrals by collisions with ions. (Further discussion on this point can be 
found in § 2.1 of Mouschovias & Ciolek 1999.) For the ion number density we assume a power-law 
behavior of the form 

- = K (ira)* • < 12 > 

where K, (~ 3 x 1CT 3 cm~ 3 ) and k (~ 1/2) are constants. In reality, the exponent k is also a function 
of density (Ciolek & Mouschovias 1998), due to the fact that ambipolar diffusion can deplete the 
abundance of dust grains in a contracting core (Ciolek & Mouschovias 1996), which alters the rate 
of capture and recombination of ions and electrons on grain surfaces. However, we ignore this effect 
in our models for the time being. 

Integrating the z-component of the force equation for the neutrals from z = to z = +Z, and 
requiring that there be hydrostatic equilibrium along field lines yields 

2 N 



B\ + B\ + 

x,Z y,Z 



R 9Z_ , r> dZ_ 

D - z dx ~r y,z dy 



p n C 2 = ^Ga 2 n + P cxt + ^ — , (13) 

where we have used the Gaussian relation for a thin sheet, g z (x,y,+Z) = —27rGa n (x,y,+Z). As 
discussed in CM93 and BM94, the first two terms on the right side of equation (13) represent the 
self-gravitational stress and external pressure acting on a sheet, respectively. The latter term - 
not included in our earlier studies (e.g., see eq. [26] of CM93) — is the total magnetic "pinching" 
term due to magnetic pressure and tension stresses that act to compress the sheet. Although this 
last term is generally smaller than the others, we retain it nevertheless in our models, since it costs 
very little computationally to include it. 

We use Poisson's equation for a very thin sheet to solve for the gravitational potential ip: 

V 2 iP(x, y, z) = 4nGa n {x, y)S(z). (14) 

Imposing the boundary condition lim^ z ^ oc Tp(x,y, z) — > 0, equation (14) can be solved by the 
method of Fourier transforms (e.g., Wyld 1976; Byron & Fuller 1992). Doing so, one finds for 

z = 0, 

F[4>(x,yM = -^G :F[a f> y)] , (15) 

K z 

1 

where T[f] is the two-dimensional Fourier transform of the function /, and k z = (k 2 . + k 2 )^ is 
a function of the planar wave numbers k x and k y . Hence, at any time t we can determine ^[tp] 
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by calculating ,F[<7 n ]. Inverting the transform then yields ip at that time, and from that we can 
calculate the gravitational field 

&fe»,0)=-^, (16a) 

%( x, 9 ,o) = -M|M. (16b) 

The magnetic field components B x z and B z can be gotten in a similar fashion. Above the 
sheet, the magnetic field can be written in the form B(x,y,z > +Z) = B'(x,y,z) + B TC {Z, where 
B' is the reduced magnetic field. The assumption that the constant pressure external medium is 
hot and tenuous implies that it is also current-free (j ext = 0). Consequently, Ampere's law gives 
V x B = in the region above the sheet. Therefore, B' = — V^, where * is a scalar magnetic 
potential, that, because of Gauss's law (V ■ B = 0), satisfies Laplace's equation, 

V 2 ^(x,y,z > +Z) = 0. (17) 

From our earlier comments on the external field we have B' —> very far from the sheet. Therefore, 

lim ■$(x,y,z) = . (18) 

z— >00 

The boundary condition on ^ at z = +Z is derived from the continuity equation (3) for the normal 
component of the magnetic field across the top surface of the sheet, which is, written in terms of 
*, 

d*(x,y,+Z) d*(x,y,+Z) dZ d*(x,y,+Z) dZ 

Wz Wx & dy ^ = -i^(x,y)-B rci ] . (19a) 

In the thin-sheet approximation, this becomes 

+ hm o d*^' +Z) = - [B z , cq (x, y) - B Ici ] . (19b) 

The solution of the magnetic Laplace equation (17) for ^, subject to the boundary conditions 
(18) and (19b), can also be performed by Fourier transform, in analogy to what is done for the 
gravitational potential ip. In the limit +Z —> 0, we find 

F[y(x,y,0)] = . (20) 

k z 

Once ^ is obtained, by inverting the transform, it follows that 

B . = lim (21a) 



dy(x,y, +Z) 
+z^o dy 



b v , z = ,k- t: t ; • (21b) 
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To close our system of equations, the evolution of the equatorial magnetic field in the plane 
of the sheet B z ^ eq (x,y) is governed by the magnetic induction equation. For the density range we 
consider in our cloud models, 10 3 cm~ 3 < n n < 10 s cm -3 , the magnetic field is effectively "frozen" 
in the ion-electron plasma (see §§ 2.2 - 2.4 of Mouschovias & Ciolek 1999 for a discussion). Hence, 
advection of magnetic flux in our system is described by 

OB d d 

-^ = -^(B z , eq v- hX )- — (B z , eq v- hy ) . (22) 



2.2. Boundary Conditions, Uniform Background State and Initial Conditions 

As described in the preceding section, model clouds are assumed to be isothermal thin planar 
sheets of infinite extent. We follow the evolution in a square Cartesian region of size L, with 
— L/2 < x < L/2, —L/2 < y < L/2. Periodic boundary conditions are used for all physical 
quantities. 

Model clouds are initially characterized by a static, uniform background state with constant 
column density a Qt o and equatorial magnetic field B rc fZ (i.e., S 2) eq,o = ^ref)- From equations 
(7a), (7b), (16a), (16b), (21a), and (21b), it is seen that all forces — thermal, gravitational, and 
magnetic — are identically equal to zero in the background state. Evolution is initiated in a cloud 
by superposing a set of random column density perturbations, Sa n (x,y) (<C cr n ,o, typically) on the 
uniform background state at time t = 0. To maintain the same local mass-to-flux ratio (T n /B Zfiq 
in the initial state as in the reference state, the magnetic field is simultaneously perturbed, with 
5B Zteq /B rc f = <5<7 n /cr n ,o- 



2.3. Dimensionless Equations and Free Parameters 

The actual system of equations we solve are dimensionless versions of the equations presented 
in § 2.1 above. We adopt the following normalizations: the velocity unit is [v] = C, the column 
density unit is [a] = cr nj o, unit of acceleration is [a] = 2irGa n fl, the time unit is [t] = C /2irGa n fl, and 
the length unit is [I] = C 2 /2irGcr n fi. From this system one can also construct a unit of magnetic 
field strength, [B] = 2irG 1 / 2 a nfi . With these normalizations, the dimensionless equations that 
govern the evolution of a model cloud are 

da n d d 

qI (ffnWn,*) = -^{vnvl^ x )- — {a n v D ,yV D , x )+a a g x -Cl a -^ + F m! , g , x , (23b) 
d . , d d i 2 ~o da a 

-^{CTnVn,y) = {^V n , y V^ x ) - — {<T n V ny ) + <J n g y - ^cfT^" + ^mag.y , (23c) 



^mag,a; — B z eq ( B _ ^ Z 



dB 



z,eq 



dx 



-12- 



mag,y 



a 



eff 



dB 



2,eq 



dt 



J i,x 



Pn 

z 

m 

9x 
On 

B x,Z 

B _ 

y,z 



ldZ 
+ 2 !te 



Bt „+Bt„ + 2B 



B 



B — Z 

y,z 



y,z 

dB 



2,eq 



dZ 



dZ 



2,eq 



dy 



idz 

+ 2~dy~ 



B\ + B 2 7 + 2B 

x,Z y,Z 



2,eq 



dx 



dZ 



dy 



B„ T — + B y Z — ] + ( £ 7 — + B„ z — ) 



<9ZY 



az 



<9x 



dy 



B„„— + B y Z — ] + [ B„ 7 — + B„ 7 — ] 



<9x 



dZ 



<9y / 



<9ZY 



dx 



y,z 



dy J 



3P cx t + <K 



2 ' 



= v n , x + 

1 

4 



(-Szjeq^ijj 



7ni,0 




<7„ 


VPn,0 


fni,0 


( Pn 


0"n 


KPn,0 



d_ 

dy 

k 



(B z ,eqVi,y) , 



F, 



mag,x j 



mag,y 



n — R — 

x > z ~dx~ + ^ z ~dy~\ 



k z 
dip 
dx ' 

dy 

-P[Pz,cq — B Ic j] 



d^> 

dx 
d^ 

dy 



(23d) 



, (23e) 

(23f) 

(23g) 
(23h) 

(23i) 

(23j) 
(23k) 
(231) 
(23m) 
(23n) 

(23o) 
(23p) 
(23q) 



where p n; o = (1 + P ex t)/4 is the dimensionless neutral mass density of the background state, and ip 
and ^ represent their values in the equatorial plane. 

The above equations are, for the most part, the Cartesian analogs of the axisymmetric thin- 
disk equations presented in CM93 (their [66a]-[66q]) and BM94 (their [34a]- [34m]). A similar set 
of non-ideal MHD equations were used to model nonaxisymmetric core formation by Indebetouw 
& Zweibel (2000), with one particular exception: they modeled clouds as infinitesimally thin, with 
Z = 0, and neglected the effect of magnetic pressure (i.e., the terms dependent on Z in the equations 
above, and also in CM93 and BM94). This means that the stabilizing effect of magnetosound modes 
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were not accounted for in the models of Indebetouw & Zweibel. Although this may be valid in some 
cases, there are, as we show below, some regions in the physically allowed parameter space for clouds 
in which magnetic pressure cannot be considered negligible compared to magnetic tension. The 
omission of magnetic pressure in these instances leads to inaccuracies in the length- and timescales 
for the onset of gravitational instability in magnetic clouds. 

Equations (23a)-(23q) have several non-dimensional parameters. P cxt = 2P ext /irGa^fi is the 
ratio of the external pressure acting on the sheet to the vertical self-gravitational stress of the 
background state. The effect of ambipolar diffusion is expressed by the dimensionless initial neutral- 
ion collision time, fn^o = 2irGa n flT n ifi/C. The limit f n i ; o — ► oo corresponds to extremely poor 
neutral-ion collisional coupling, so that the ions and magnetic field have no effect on the neutrals. 
In the opposite limit, f n i ; o = 0, the neutrals are perfectly coupled to the ions, due to frequent 
collisions, and the magnetic field will be essentially frozen in the neutral matter. The parameter k 
(=1/2, typically) is the exponent in the power-law expression (12) used to calculate the ion density 
as a function of the neutral density. As we noted earlier, this constant power-law assumption 
is only an approximation, since ambipolar diffusion has been shown to make k a function of n n 
in models with a more realistic ion chemistry network (Ciolek & Mouschovias 1998). Finally, 
B Te { = B re f/27rG 1 / 2 a Q) o is the dimensionless magnetic field strength of the background state. For 
the units we have chosen for the column density and the magnetic field, the dimensionless mass- 
to-magnetic-flux ratio of the background state is 

W = 2*0*%*=* . (24) 

-Dref -D rc f 

This also happens to be the mass-to-flux ratio in units of the critical value for gravitational collapse, 
l 

l/27r£r2 (Nakano & Nakamura 1978; see, also, § 3.1 below). Models with < 1 (B TC { > 1) are 
subcritical clouds, while those with fj,o > 1 (-B re f < 1) are supercritical. 

Normally, we set a ni o by specifying the temperature T and the density n Di o of the background 
state. Using equation (13), and the fact that C = 0.188(T/10 K) 1 / 2 km s _1 for m n = 2.33 a.m.u., 

3.63x10-^ n n , U/ T \3 _ 2 
Vn,0 = ; r ma „„- 3 77TF7 § cm ■ ( 25 ) 



(1 + P-)* Vl ° 3cm ~ 3; Vl ° K/ 

Both C and cr n; o are normalizing units that we have adopted for our cloud models. The units for 
length and time have the scalings 

and 

[t] = 3.98x10^(^^^(1 + ^)' y, (27) 
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The unit of mass is then 

3 I 



,M,^,ax 1 o-( I ^) 5 (^) 5 ( 1 + P„^M a . (*) 



It follows from equations (11), (12), and (25) that 



0.241 /3 x 10- 3 cm' 3 \ /10 5 — " 3N * ° 



(1 + P«t) 5 

From equations (24) and (25), a model's reference magnetic field is given by 

5.89 x 10~ 6 / n T1 n \ 5 / T \ 5 



5 



ref 



^(l + ^ext) 2 



2.4. Numerical Method of Solution 



The system of dimensionless partial differential equations presented in the preceding section is 
solved by the method of lines (e.g., Schiesser 1991). This was used in the axisymmetric mod- 
els of Morton, Mouschovias, & Ciolek (1994), Ciolek & Mouschovias (1994, 1995), and Basu 
& Mouschovias (1994, 1995a, b), and the nonaxisymmetric models presented in BC04. In this 
method, spatial derivatives within the PDEs are approximated by finite differences. The square 
computational domain of size L x L is discretized by dividing the region into iV 2 uniform cells of 
size L/N x L/N. L is typically chosen to be a factor ~ 4 greater than the wavelength of max- 
imum gravitational instability AT, m (see § 3.1 below). Three-point centered differences are used 
to approximate gradients. Advection of mass and magnetic flux is performed with the monotonic 
upwind scheme of van Leer (1979). The spatial discretization converts the system of PDEs to a 
system of coupled ordinary differential equations of the form dy/dt = F(y,t), where dy/dt, y, 
and F are all vectors of dimension VN 2 (V is the number of dependent variables). An implicit 
Adams-Bashforth-Moulton method is then used to time-integrate the resulting system of ODEs. 

At each time step, numerical solution of the Fourier transform and also the inverse transform of 
quantities in the equations for the gravitational and magnetic potentials ([231] and [23o]) is carried 
out with standard two-dimensional Fast Fourier Transform (FFT) techniques (e.g., Press et al. 
1996). 



3. Stability 

We now consider the response of model 
bations. Such an analysis will illuminate the 



of Cloud Models 

clouds to small-amplitude disturbances or pertur- 
basic physics of gravitational instability in weakly 
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ionized, sheet-like magnetic clouds, and will also provide the basis for understanding the length 
and time scales over which instabilities will develop in fully nonlinear calculations. 

The linear stability of isothermal self-gravitating equilibrium layers with frozen-in magnetic 
fields was examined early on by Nakano & Nakamura (1978). They determined the critical mass- 
to-flux ratio for such objects, which is discussed further below. The gravitational stability of 
weakly-ionized and thin (but with finite thickness) axisymmetric magnetic disks was investigated 
by Morton (1991) and Morton & Mouschovias (1991, unpublished). Their analysis is very similar to 
the one that we present in this paper, and we reprise many of their results in the following section. 
A later, independent linear study of partially- ionized magnetic disks/sheets was also presented 
by Zweibel (1998). Similar to that which was done in Indebetouw & Zweibel (2000), Zweibel 
(1998) studied clouds that were infinitesimally thin [Z = 0), and concentrated on those that are 
magnetically subcritical. 



3.1. Linearization and Analysis 

As is frequently done, the unperturbed zero-order (background) state of a model is assumed 
to be uniform and static. The dimensionless equations presented in § 2.3 above are linearized to 
first-order by writing, for any physical quantity, 

f(x,y,t) = f + Sf 

= f + Sf li e i ^ x+k y y - ojt \ (31) 

where /o refers to the unperturbed state, Sf is the perturbation and <5/ a is its amplitude, k x and k y 
are again the x and y wave numbers, and to is the complex angular frequency. The perturbations 
are small, with \5f a \ <C fo- With regard to velocities, which have vq = because of the static 
assumption, it is understood that \5v\ <C characteristic signal speeds of the multifluid system (i.e., 
sound speed, Alfven speed, etc.). 

For the assumed type of perturbation (31), we can set d/dt — > —iuj, d/dx — > ik x , and d/dy — > 
iky. The dimensionless equations for a model cloud then become, collecting and retaining terms to 
first-order, 

lo 8a n = k x 8v QtX + k y 5v nw , (32a) 
lo 5v ntX = y (c^ Sfi k z - I) 5a n + (1 + k z Z ) 5B Zfiq , (32b) 

uSv Q>y = ^{c 2 cSfl k z -l)5a n + B rci ^{l + k z Z )5B ZfiCl , (32c) 
w SB Zteq = B rei k x 5v ntX + B iei k y Sv n: y - if nifi B? ci k z (1 + k z Z ) 5B Zfiq . (32d) 

From equation (23f), the reference state effective isothermal speed of sound is 

l 

(l + 3P cxt ) 1 

Ceff '° = — v~ ~r> — • ( 33 ) 

-L + -Text 
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In deriving the linearized system (32a)-(32d), we have used equations (231)-(23n) to relate gravi- 
tational field perturbations to column density perturbations, (23o)-(23q) to relate planar magnetic 
field perturbations to those of the equatorial vertical magnetic field, and (23h) and (23i) to sub- 
stitute the perturbed ion velocity components with those of the neutrals. The wavenumber k z has 
the same meaning as previously defined. A mode is unstable if the imaginary part of the complex 
frequency u> l > 0. The growth timescale of the instability is r g = l/wj. 

The fundamental physics of the linear system is readily discerned from the various terms in 
equations (32a)- (32d): thermal-pressure and self-gravitational forces are proportional to perturba- 
tions in the column density in equations (32b) and (32c), and magnetic tension and pressure forces 
are proportional to perturbations in the equatorial magnetic field. The drift or diffusion of magnetic 
field and plasma with respect to the neutrals is represented by the term in the magnetic induction 
equation (32d) that contains f^o; comparison with the linearized mass continuity equation (32a) 
shows that, in the limit f n i ; o = 0, collisional coupling of the neutrals and ions (and therefore, the 
magnetic field) is instantaneous and perfect, and they all move together as a single fluid. In the 
opposite extreme, f^o 3> 1, neutral-ion collisions are infrequent, and the neutral and ion-magnetic 
field fluids are increasingly decoupled and move independently of one another. 

The solution of the full dispersion relation for the gravitationally unstable mode of the lin- 
earized equations (32a) - (32d) are presented in Figure 1 for (a) fn^o = 0, (b) fn^o = 0.04, (c) 
7ni,o = 0.1, (d) fn^o = 0.2, (e) fni,o = 1, and (f) fni,o = 10. The external pressure parameter for 
these models is P cx t = 0.1, which sets the dimensionless effective sound speed C e fr,o = 1-04. (Note 
that CeS,o = 1 at P cxt = and 1. Ceff.o is maximal at P cxt = 1/3, and is equal to 1.061.) Displayed 
in each panel of Figure 1 is the growth time r g as a function of the wavelength A (= 2ir/k z ), for 
various values of jjlq = l/B re {. The separately labeled curves show the result for //q = 0.5, 0.8, 1, 
1.1, 2, and 10, respectively. We are able to use A as the independent variable because the char- 
acteristic polynomial for our eigensystem is found to be only a function of k z = (A; 2 . + ky)% . This 
means that to this order of approximation, all perturbations are independent of the planar angle 
of propagation 9 [= taxL~ 1 (k y /k x )]. 

Understanding the data presented in Figure la - / is aided by examining the results for the 
instability growth timescale in the following two limits. 



In this limit, the neutrals, ions, and magnetic field respond as a single combined fluid, with 
the magnetic flux frozen in the neutrals, and the resulting dispersion relation is found to be 



The gravitationally unstable mode corresponds to one of the roots of uj 2 < 0, and occurs for 



Limit 1: Flux- freezing, f n i 5 o = 0. 




(34) 
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B rc f < 1, or, equivalently, /xq > 1. The growth timescale for this mode can be written as 



A 



27r(l-4 2 cf )(A-A M s) 



l ' 

2 



(35) 



for A > AmS) where 



Ams = 2-7T 



C eS,0 



B ?cf Z 



1 



rcf 



(36) 



is the critical or threshold wavelength for instability. The minimum growth timescale (= maximum 
growth rate) occurs at AMS,m = 2Ams- The flux- freezing growth timescale (35) for our various 
model clouds is plotted as a function of A in Figure la - /, shown as open circles. 

We have given to the critical wavelength Ams the subscript 'MS' because it is the maximum 
lengthscale that can be supported by both magnetic and thermal pressure effects, and therefore is 
related to magnetosound modes. This can be readily seen by noting that the dimensionless Alfven 
speed in our model clouds is given by 



Va,o 



B rc { 
\/2pn,0 



(37) 



(see eq. [23k]) and the fact that the dimensionless column density is unity for an unperturbed 
cloud. It follows then that the isothermal magnetosound speed (since we consider only isothermal 
perturbations in our models) in the adopted set of units is 



^MS,0 = (C'eff.O + ^A,o) 2 

1 

2 (1 + Pext) 



C e fT,0 



1 + 



Mo (1 + 3Pext) 



In the brackets of the last equality we have used eq. (33) and the relation 

2 



(1 + P ( 



(38) 



(39) 



cxt J 



which follows from the linearization of equations (23j) and (23k), to eliminate C^ s and Zq in terms 
of P e xt- Examination of equation (36) reveals the presence of the magnetosound speed (38) in this 
expression, thus identifying the combined action of thermal and magnetic pressure in the support 
of a cloud against self-gravity. It also shows the importance of magnetic pressure in setting the 
instability timescale and the wavelength of maximum growth rate. For clouds that are close to 
being critical (/^o ~ 1), neglecting this effect (e.g., Zweibel 1998; Indebetouw & Zweibel 2000) can 
significantly underestimate r g and Ams,hi- 
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When there is negligible magnetic support, that is, when £? re f — > (/xo — > °°) it follows that 
Vms,o ->■ C'cff.o, and A M s -> A T , where 

A T ee 27 r(7 e 2 ffi0 = tt (^2=^ (40) 

\ 1 + "ext / 

is the critical thermal lengthscale. For this situation, the growth timescale for the unstable mode 
is still given by equation (35), with B ve { = 0, and At replacing Ams- The maximum growth rate 
for the unstable mode in this circumstance then occurs at AT, m = 2 At- 

Limit 2: Stationary magnetic field lines, uSB z>eq = 0. 

For this situation, which would be relevant to models with effective ambipolar diffusion, and 
therefore, relatively weak coupling of the neutrals to the ions and magnetic field, the resulting 
dispersion relation is 

u, 2 + ^-u - ((7 C V 2 - k z ) = 0. (41) 

T n i,0 V J 

From this relation, one finds that an unstable mode exists for A > At, and has a growth timescale 

r g = ^ — • (42) 

2 



A 2 + 8vrf 2 (A-A T ) 



- A 



The minimum growth time for this mode is at the wavelength AT, m j defined above. The growth 
time (42) is also displayed as a function of wavelength in Figure la - f (crosses). 

As f n i 5 o — > oo, r g — > X/y / 2ir(X — At), which is identical to the result from equation (35) when 
Bid = 0. This is because in both of these circumstances the magnetic field does not affect the 
neutrals: fni,o — ► co corresponds to when there is no collisional coupling between the neutrals 
and ions (and, hence, the magnetic field). The ions are completely "invisible" to the neutrals 
in this situation, and there is no transmission of magnetic force to them via neutral-ion collisions. 
Similarly, when B TC { = 0, there is no magnetic field, and therefore magnetic forces do not contribute 
to the support or dynamics of a model cloud in that case. 

With the results of the two limiting cases 1 and 2 described above in hand, the underlying 
physics of the gravitationally unstable modes presented in Figure la - / is made more transparent. 
For instance, the models with the magnetic field frozen into the neutrals (f^o = 0) displayed in 
Figure la are seen to be in exact agreement with the growth timescale predicted by equation (35). 
Particularly noteworthy is the dependence on the parameter (= l/i? re f) for these models: there 
is no unstable, gravitationally collapsing mode for [1q < 1. Moreover, models with (1q ~ 1 have 
minimum growth timescales at wavelengths AMS,m that are much larger than those for models with 
greater values of fiQ. As an example, we note that Ams.hi for the model with fio = 1.1 in Figure la is 
an order of magnitude greater than that in the model with = 2. The actual value of the growth 
timescale for the fio = 1.1 model is seen to be about an order of magnitude greater than that for 
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the hq = 2 model as well. This is a consequence of there still being near-equality of gravitational 
and magnetic forces when is very close to unity. 

When there is imperfect neutral- ion coupling (f^o > 0), gravitational instability is possible 
even for model clouds with fiQ < 1. This can be seen in panels b - f of Figure 1. The timescale for 
the instability is typically greater than that for the supercritical flux- frozen models (Fig. la), but 
it is finite. For these models, a Jeans-like growth of density perturbations occurs with the collapse 
moderated by the retarding collisional forces exerted on the neutrals as they diffuse through the 
plasma and field. A result of this type was first noted by Langer (1978). Models that are more 
subcritical are better approximated by the stationary field limit described above. This can be seen 
by the excellent agreement of the limiting growth timescale given by equation (42) — and also 
the wavelength of maximum growth (Ax,m) — with the dispersion curve for the fio = 0.5 model 
displayed in Figure lb - f. As /iq approaches unity, there is a transition of the growth timescale 
behavior from the stationary field/ambipolar diffusion limit (eq. [42]) and toward the flux- freezing 
limit (eq. [35]). For hq close to unity, the instability proceeds through a hybrid mode in which 
both ion-neutral drift and gravitational contraction with field-line dragging are active. As discussed 
above, the two limiting approximations approach one another for models with fni o 3> 1 and ^ > 1, 
since either limit describes clouds with effectively no magnetic support. This accounts for the near 
equality of all the model curves seen in Figure If. 

Figure 2 presents the wavelength A giIn , which is defined as the wavelength A that has the 
minimum growth time, as a function of /j,q f° r models with the same values of f n i ; o and P cxt as 
in Figure la - f. For comparison, the dashed line in Figure 2 is the constant value of Ax.m for 
those particular models. Consistent with our discussion above, we note that there is singular and 
limiting behavior of this lengthscale for the model cloud with fn^o = and fio = 1. The value of 
Ag im is also seen to be especially sensitive to the value of Ho for near-critical clouds with f^o > 0, 
exhibiting a sharp, resonant-like peak in the region /io ~ 1- Table 1 lists the ratio Ag )m /Ax,m at 
the peak for each of the non- flux- frozen models in Figure 2. The peak ratio is largest for models 
with f n i 5 o < 1, while, for models with f^o 3> 1, magnetic field effects are barely transmitted to the 
neutrals (because neutral-ion collisions are very rare), and this ratio instead approaches unity. 

Based on these results, we expect then, that clouds that are marginally or slightly supercritical 
will form cores that evolve more slowly, and have size scales (radii, and core spacing) significantly 
greater than that which will occur in clouds that are more highly supercritical. In addition, nearly 
critical clouds will also have size scales that are markedly larger than those in distinctly subcritical 
clouds. As we shall see in a following nonlinear study, this specific dependence of the gravitational 
instability on the initial mass-to-flux ratio leads to notable differences in the physical characteristics 
of collapsing cores in magnetic interstellar clouds that should be easily discerned by observations. 

Figure 3a - f shows the growth timescales for the same model clouds as in Figure 1, but with 
the dimensionless external pressure P cxt = 10. In the limit of large P ext , equation (23j) shows that 
the cloud is pressure confined, with Z oc a D (eq. [23k]). As a result, a local peak in a D is a peak in Z. 
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For this situation then, the top surface of the cloud looks like a dome and the external pressure force 
has horizontal components that point inward to the dome's peak, acting in the same directions as 
the gravitational field components g x and g y . This further enables the gravitational clumping, and 
decreases the growth timescale, appearing as a reduction of C e g t o in our equations. Because of this, 
the dimensionless sound speed C^o is decreased by a factor of 2.05 from its value in the models 
with P ex t = 0.1. This also reduces the growth timescale r g and associated lengthscales Ams and At, 
since they are also functions of C e ff,o, an d therefore, of P cx t- The net effect is to shift and reduce 
the instability timescale and critical wavelengths by a corresponding factor from those seen in the 
models with lower external pressure. This means that clouds in regions with much greater external 
pressure (perhaps due to being embedded in a massive cloud complex or fragment, or adjoining 
hotter environments such as an HII region and/or shocked gas) will have characteristic sizes that 
are smaller when compared to clouds surrounded by lesser external pressures. However, the general 
behavior and physics of the models as a function of Tnio and fj,Q, including the applicability of 
the limiting analytical approximations discussed above, is similar to that seen in the previously 
described models with smaller P e xt- 

The wavelength with the minimum growth time A giI11 is shown as a function of jiQ in Figure 4 
for models with the same values of fni o and P e xt as in Figure 3. It is again the case that, for clouds 
that are almost critical and with fni,o % 1, the wavelength that has the most rapid gravitational 
response is significantly greater than in models that are farther away from fiQ = 1. The peak ratio 
values of A giI11 to AT, m for the models with fn^o > are also listed in Table 1. As before, this ratio 
is much greater for the models with f^o < 1. 

3.2. Comparison to numerical simulations 

We now compare the results of the dispersion analysis of the preceding section to simulations of 
the early evolution of model clouds governed by the full system of non-ideal magnetohydrodynamic 
equations (23a) - (23q). This system of equations is solved by a numerical code using the techniques 
described in § 2.4. The code was previously used to generate results for two model clouds whose 
nonlinear evolution was presented by BC04. During the early evolution of a model cloud, physical 
variables do not change very much from their initial state values. Hence, during this period the 
physical evolution should be close to that determined by the linear analysis above. Comparing the 
results of the early-time evolution of the full simulations to those of the dispersion analysis allows 
us to reinforce our understanding of the underlying physics governing the early evolution of a cloud 
— the precursor to the later nonlinear phase of evolution — as well as provide a useful benchmark 
and theoretical standard to establish the overall accuracy of our numerical code. 

Because we are concerned with the linear stages of gravitational instability, we will confine our 
focus of the numerical simulations in this study to the evolution of the column density in model 
clouds. The detailed time-dependent behavior of other quantities, such as the velocity fields of 
the ions and neutrals, and the evolution and redistribution of mass in magnetic flux tubes will be 
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described in a following paper. 



3.2.1. Monochromatic perturbations 

We follow the evolution of clouds initially given perturbations with a single wavelength A. The 
evolution of a model cloud is initiated by imposing at time t = an initial (normalized) column 
density profile of the form 

a n (x,y,0) = I + 5o- nA cos(2irx/\), (43) 

that is, a uniform background state with a perturbation with amplitude 5a Iha . In using this 
particular perturbation, we make use of the fact that, as mentioned in § 3.1 above, the dispersion 
analysis indicates that the effects of linear disturbances are independent of the angle of propagation 
8. Since it really doesn't matter for our purposes which direction of propagation we choose, we 
take 6 = (parallel to the x-axis), which means that k y = 0, and therefore k x = k z = 2tt/X. The 
initial velocity and magnetic field are also perturbed in a way that is consistent with the system of 
equations (32a) - (32d) for the column density perturbation specified by equation (43). Solving for 
the initial perturbations Sv n>x , 5B Zteq , and 5v^ y in terms of the given 8a n , k x , and r g = i/co (which 
is a function of A) from this set of equations yields 

Sv n x(x,y,Q) = -— — <5o- nia sin(27rx/A), (44) 
2vrr g 

\B f 

SB Zteci (x,y,0) = — — z rcf 5a nta cos(2ttx/X) , (45) 

A + 27TTgT nii o-B i r cf (1 + 2-kZq/X) 

and 8v n! y(x, y, 0) = 0. By specifying the relation between the perturbed physical variables in this 
way, we are basically selecting the eigenvector of the perturbation at the wavelength A. Hence, the 
monochromatic perturbation excites a single eigenmode of a model cloud at t = 0. 

Not surprisingly, when we initiate the time evolution in this fashion, the subsequent evolution 
of a model cloud is simply the continued growth of the excited lone eigenmode. The left panel 
of Figure 5 displays the growth of the column density maxima for four different cloud models as 
a function of time. Each has f n i ; o = 0.2 and P ext = 0.1. For all of these models, 5<r n ,a = 0.02 
and A = 4ir. The calculations were performed on an equally spaced mesh of 32 x 32 cells on 
a square computational region of size L = 4n in each direction. Shown in that Figure (solid 
curves) is the evolution of the maximum reduced column density 6o~ n (t) = o~ n (t) — 1 for models 
with fio = 0.5, 1, 2, and 10. Also shown (dashed curves) is the growth of the column density 
perturbations that would occur as predicted by the linear relation for an unstable mode at the 
same location, 6a n (t) = 5a D ^exp(t/r g ). The values of r g used to plot the curves are those derived 
from the earlier dispersion analysis for A = 47r, f^o = 0.2 and P ex t = 0.1. From the data in Figure 
Id we find that r g = (20.3, 16.0, 4.34, 2.12) for fi = (0.5, 1.0, 2, 10). Comparing the numerical 
simulation results to the theoretical values predicted by the linear analysis indicates that they are 
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in excellent agreement and overlap during much of the early evolution of each model. In fact, 
significant deviation (> 1%) between the simulation results and the linear predictions does not 
begin to occur until 5a a has grown to ~ 0.2, which is well beyond what is generally considered 
the regime of linear growth. At later times nonlinear effects are evident, and the column density 
grows more rapidly and significantly exceeds that predicted by the linear analysis by the end of the 
period shown in Figure 5, when 8a n = 1. 

The growth of the column density perturbations in the numerical simulations at early times 
is entirely consistent with our discussion of the physical processes acting in the dispersion analysis 
in § 3.1. For instance, the linear analysis predicts that the onset of gravitational instability in the 
subcritical /io = 0.5 model is due to the action of ambipolar diffusion (see Fig. Id). Moreover, 
using the stationary field limit approximation (eq. [42]) to calculate r g for this model (At = 6.75 
for these parameters) gives r g = 21.8, which differs from the full dispersion analysis result quoted 
above by only 7%. The critical model = 1 and supercritical model fio = 2 are also heavily 
influenced by ambipolar diffusion, and they would not be able to collapse in its absence, because 
the excited wavelength for both of these models is in the region At < A < Ams- Thus, these 
two models lie in the transition region between diffusion-regulated and flux-frozen collapse for this 
particular wavelength, as noted in our discussion of these models in regard to Figure Id. Finally, 
the instability of the highly supercritical model = 10 occurs essentially with freezing of magnetic 
field lines in the neutral matter. For this model, Ams = 6.93 which is < 4ir. Using the flux-frozen 
approximate expression for the growth time given by equation (35) yields r g = 2.12 for this model, 
which is exactly the same as the full dispersion relation solution, and is also in very good agreement 
with the model simulation's early-time linear growth. 

The right panel of Figure 5 shows the growth of the column density maxima in another four 
model clouds, with the same wavelength and parameters except P e xt = 10. The development of 
these higher- P ext models is again seen to be in excellent agreement with the dispersion analysis for 
the linear stage of evolution. They are also in accordance with the values of r g predicted by the 
analytic relations (42) and (35) for the subcritical and supercritical models, respectively. 

We conclude then, that our numerical code accurately reproduces the physics of gravitational 
instability in planar magnetic model clouds. 

3.2.2. White noise perturbations: the wavelength of minimum growth time A gim and its 

dependence on hq 

In this subsection we consider cloud models that are given a spectrum of random, small- 
amplitude perturbations in the physical variables at t = 0. The spectrum we use is white noise, i.e, 
flat, so that there is no preferred wavelength selected by the perturbations. However, we introduce 
damping so that wavelengths equal to twice the mesh spacing and smaller are negligible. The 
root-mean-square amplitude of the fluctuations in the initial state of a model cloud is 3% of the 
uniform background. The size of the computational region in each direction is L = IQtt, and the 
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number of cells in each direction is N = 128. 

The white noise spectrum will excite an ensemble of eigenmodes with a wide range of wave- 
lengths. The dominant evolutionary modes that will emerge from this ensemble of fluctuations and 
govern the subsequent evolution will be those grouped about the one with the minimum growth 
time. A cloud will then develop features with lengthscales comparable to the wavelength with the 
minimum growth time, A gim . That this intuitive notion is valid can be seen in Figure 6, which 
shows contour maps of the column density a D (x, y) in two model evolution simulations at the time 
for which the value of the maximum column density has grown to cr njUiax = 2 (twice the initial 
background value) in each model. Both models again have f^o = 0.2 and P ext = 0.1. The model 
displayed in the left panel is initially subcritical with fio = 0.5, and the one in the right panel started 
as a critical cloud with fiQ = 1. The density contours cover the range 0.8 to 1.6, in increments of 0.2. 
By the time that the contour snapshots are taken, nonlinear effects have set in. This is evidenced 
by column densities that are relatively large compared to that of the initial background state, and 
the beginning of nonaxisymmetric gravitational fragmentation (which will form protostellar cores) 
due to the interaction of various eigenmodes. Despite this, both cloud models retain a significant 
imprint from the predecessor linear stage of evolution: there is a nearly uniform separation between 
core fragments — defining here a core as the enclosed higher column density regions with a n > 1.2 
- and the fragmentation scale for these nascent cores is close to A gim in each model. For the model 
with fiQ = 0.5, Figure 2 indicates that A giI11 = 14, while for the cloud with /io = 1> A giI11 = 24. Exam- 
ination of Figure 6 reveals that these scales correspond to the mean distance between the different 
a n = 1.2 contours for these models. Additionally, we note that the behavior of the fragmentation 
scale in these models as a function of /io is also consistent with the linear analysis. Notably, the 
average spacing between cores in the critical model evolution simulation is larger than that in the 
subcritical model, in line with the predicted dependence of A gjm on /iq (see Fig. 2). This results in 
there being fewer cores within the same size region. In fact, according to Table 1 (see, also, Fig. 
Id), the maximal A giin (= 52.8) for these values of fni,o and P ext occurs at fio = 1.13. 

Further confirmation that A giin sets the characteristic lengthscales in clouds can be seen in 
Figure 7, which presents the column density plots of two more model clouds, with the same values 
of f n i 5 o and P ext as in Figure 6. The model in the left panel of Figure 7 is supercritical and has 
Ho = 2. The linear analysis predicts A g)m = 23 for this model, which is indeed close to the the 
distance between core-bounding contours in that Figure. There are also relatively fewer cores that 
have formed, similar to that seen in the critical model. The similarity in the structural (lengthscale) 
properties of the = 1 and hq = 2 models is not surprising, as they have nearly equal values for 
A g)in . However, the dynamical properties of the = 1 and /xo = 2 models are quite different in 
the nonlinear regime, as shown by BC04. Specifically, the maximum velocities in the supercritical 
models become supersonic on scales that are well within the resolving power of modern observations. 
We defer the detailed discussion of the dynamical and related nonlinear evolution of cloud models 
to a study that will follow this paper, and to that already presented in BC04. Finally, the right 
panel of Figure 7 shows the column density contours for a highly supercritical model with = 10. 
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For this model, A g)II1 = 13. This is also in agreement with the mean distance between cores in the 
contour plot. The lengthscales and number of cores in the highly supercritical model are akin to 
that seen in the subcritical fiQ = 0.5 model. This is no coincidence, since these two models have 
similar values of A giII1 . Although a somewhat counterintuitive result, this is because there is also 
a Jeans-like growth of density perturbations in the subcritical models, but it occurs on the much 
longer ambipolar diffusion timescale as neutrals drift past near-stationary plasma and magnetic 
field. 

There are also some differences in the models in Figures 6 and 7 that become more pronounced 
with increasing initial mass-to-flux ratio: the cores become more nonaxisymmetric with greater 
values of /j,q, showing enhanced elongation along a single axis. This second-order effect — an 
amplification of nonaxisymmetry due to the nonlinear interaction of certain eigenmodes for the 
more supercritical models — will be a topic of investigation in a following paper. And, as mentioned 
above, the dynamical properties of the model clouds, such as their velocity fields, are also markedly 
different, characteristically having greater infall velocities over larger scales with increasing jiq. 
This too will be studied further in a future paper. Despite these differences, we conclude from our 
representative evolution simulations that the gravitational fragmentation scale of a model cloud 
with an initial spectrum of perturbations is effectively the value of A gjm for that particular model. 
The fragmentation scale has a dependence on hq in agreement with the result of the linear analysis, 
especially with the prediction that near the critical value there is a dramatic increase in A gjm . 
Because of this resonant-like behavior, clouds that are close to critical will have significantly larger 
size scales than in models that are more highly supercritical or more highly subcritical. 

4. Summary 

We have presented the formulation of physical models that we will use to study the nonax- 
isymmetric formation and self-gravitational collapse of protostellar cores in magnetic interstellar 
molecular clouds. Model clouds are partially ionized isothermal thin planar sheets with finite 
half-thickness Z(x,y). 

The system of equations that are used to govern the time evolution of model clouds contain four 
fundamental parameters, all defined in § 2.3. The first is the dimensionless background magnetic 
field strength of a cloud, B re f, the inverse of this parameter is the initial mass-to- magnetic flux 
ratio in units of the critical value for gravitational collapse, /j,q. Clouds with B TC { > 1 (/xo < 1) 
are subcritical and are initially magnetically supported, and cannot collapse in the absence of 
ambipolar diffusion. Models with B TC { < 1 (/xo > 1) are supercritical from the outset and unable to 
support themselves against their own self-gravity. The next parameter is the normalized neutral-ion 
collisional momentum-exchange time f^o- It is a measure of the efficiency of ambipolar diffusion in 
model clouds: f^o — ► corresponds to very effective momentum transfer and collisional coupling 
of neutrals and ions, resulting in freezing of magnetic flux in the neutral matter. In the opposite 
limit of fn^o S> 1, collisions are so infrequent that the neutrals are essentially decoupled from the 
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plasma and magnetic field, and magnetic forces contribute negligibly to the support and evolution 
of a model cloud. P cxt is the ratio of the external pressure to the vertical self-gravitational stress 
in the initial uniform reference state of a model cloud. The final parameter is the exponent k in 
the power-law expression used to calculate the ion density as a function of the neutral density. 

We also investigated the linear stability of model clouds, and how the growth times r g and 
critical wavelengths of the gravitationally unstable modes depend on /io, fn^o, and -P ex t- Analytic 
expressions for r g and the critical wavelengths were derived for the limits of weak and strong 
ambipolar diffusion. These expressions (eqs. [35] and [42]) agreed well with the full dispersion 
results in these limits. 

Models with frozen-in magnetic flux (fn^o = 0) are gravitationally unstable when they are 
supercritical and the wavelength of a perturbation exceeds the magnetosound critical wavelength 
Ams- For finite thickness clouds, magnetic pressure contributes substantially to setting the value 
of Ams an d r g when /io ~ 1. Ambipolar diffusion (f n i,o > 0) allows clouds to be unstable even 
when subcritical, so long as the perturbation's wavelength is greater than the thermal critical 
wavelength At (At < Ams> scc c q s - [36] and [40]). The instability modes of clouds with P ex t > 1 
behave qualitatively the same way as for those with P cxt < 1, except that quantitatively their 
critical wavelengths and growth times are reduced, due to the effective retarding sound speed being 
decreased at higher external pressures. The dispersion analysis also revealed how the wavelength 
with the minimum growth time, A giin , behaves as a function of /io- For models with T n [ ? o ^ 0, A g?m 
has a resonance at a value of /io that is usually just slightly greater than the critical value /io = 1- 
Because of this resonance, a model cloud with /io ~ 1 will have a A gim that is significantly greater 
when compared to that of models that are much more sub- or supercritical. 

In addition to providing basic insight to the physics of gravitational instability in partially 
ionized media, the linear analysis was also used as a theoretical standard to test the accuracy of 
our numerical method of solution of the full set of nonlinear equations (eqs. [23a] - [23q]) that 
govern the time evolution of a model cloud. The early-time evolution of several cloud simulations 
given an initial monochromatic perturbation of wavelength A was compared to the predicted value 
of the growth timescale from the dispersion analysis. The temporal growth of the column density 
maxima in each model simulation was in excellent agreement with the theoretically predicted values 
of r g , thus establishing the accuracy of our non- ideal magnetohydro dynamic computational code. 

Finally, we presented a suite of model simulations that had their evolution initiated by fluctu- 
ations with a spectrum of wavelengths, following their evolution a little beyond the linear growth 
phase. The characteristic fragmentation scale that developed in these models tends to correspond 
to the wavelength with the minimum growth time A giin of the initial state. The resonance behavior 
of Xg t m for clouds with the parameter /io near the critical value — as predicted by the linear analysis 
- was also in evidence in the numerical simulations. Those having /io ~ 1 had a much greater 
mean spacing (« A gim ) between cores than in models with /io <C 1 and /io S> 1. Because of this, the 
total number of cores that developed in the models near the resonance was less than in those that 
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were farther away. This sensitive dependence of A g)m about (io ~ 1 has important implications for 
core and star formation, since observations currently indicate that mass-to-flux ratios in molecular 
clouds generally lie in the range 0.5 < no < 2 (Crutcher 2004). 

In forthcoming studies we will explore further the formation and nonaxisymmetric collapse of 
protostellar cores as a function of the fundamental model parameters, building on the formulation 
and analysis presented in this paper. The properties of dynamically infalling cores such as spatial 
density and velocity maps, core shapes, magnetic field strengths, and other quantities of interest, 
will be investigated in detail. In doing so, we will be providing a physically consistent model with 
testable, quantitative predictions that may be used to interpret and perhaps guide observations of 
star formation in magnetic interstellar molecular clouds. 
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neering Research Council of Canada. Helpful comments from an anonymous referee are gratefully 
acknowledged. 



-27- 



Table 1 

Peak ratio of wavelength with minimum growth time (A gjm ) to wavelength of 

MAXIMUM GRAVITATIONAL INSTABILITY (A Tjm ) 



P cx t =0.1: 



Pcxt = 10: 
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Fig. 1. — Instability growth times for various model clouds as a function of wavelength, for models 
with (a) f n i,o = 0, (6) f n i,o = 0.04, (c) f n i,o = 0.1, (d) f n i,o = 0.2, (e) f n i,o = 1, and (J) f n i,o = 10, 
respectively. In each panel are shown the timescale curves for models with mass-to-flux ratios 
Ho = 0.5, 0.8, 1, 1.1, 2, and 10 (labeled). Each model has P ex t = 0.1. Also shown are the results 
of the approximate analytical solutions in the limit of flux-freezing (open circles), and stationary 
magnetic field lines with ambipolar diffusion (crosses), given by eqs. (35) and (42), respectively. 
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Fig. 2. — Wavelength with the minimum growth time (= maximum growth rate) as a function 
of initial mass-to-flux ratio. The displayed curves are for models with f n i,o = 0, 0.04, 0.1, 0.2, 1, 
and 10, respectively. All models have P cxt = 0.1. The dashed curve is the value of AT, m for these 
models. 
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Fig. 5. — Time evolution of column density maxima for numerical simulations of various cloud 
models given an initial monochromatic perturbation of wavelength A = 4n. The amplitude of 
the initial perturbation of the column density is 5a nA = 0.02. Each model has f n i ; o = 0.2, and 
is labeled by the value of its mass-to-flux ratio fiQ, which is 0.5, 1, 2, and 10, respectively. Left 
panel: models with P ext = 0.1. Right: models with P cxt = 10. The solid lines show the peak 
density in each model simulation. The dashed lines are the theoretical values given by the relation 
S(T n (t) = Sa n>a exp(i/r g ), where r g is the growth time for each model predicted from the linear 
analysis. 
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Fig. 6. — Column density contours of two model clouds that initially had random small- amplitude 
fluctuations in physical variables superposed on a uniform background state. The contours arc 
overlaid on a grayscale image of the logarithm of the column density. At the time shown, the peak 
column density in each model is twice the value of the background state (<T n ,max = 2). Each model 
has f n i 5 o = 0.2 and P ex t = 0.1. The left model has /xo = 0.5, and the right one has the critical value 
= 1. The density contours range from 0.8 to 1.6, in steps of 0.2. 
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Fig. 7. — Same as in Fig. 6, but for supercritical models with /j,q = 2 (left) and no = 10 (right). 



